clear
#delimit ;
set more off;

use pw_coalitions.dta;

gen logshift=ln(1+maxshift);
gen logloser=ln(1+losercinctv);
gen logint=logshift*logloser;

fp <pceyr>: probit failure logshift logloser logint lognum loglength <pceyr>, cluster(warnum);

preserve;
set seed 304876538;

drawnorm SN_b1-SN_b8, n(10000) means(e(b)) cov(e(V)) clear;
postutil clear;
postfile mypost2 prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi
using pwsims2.dta, replace;

noisily display "let's do this!";

local a=0 ;
while `a' <= 0.282714 { ;

{ ;

scalar h_logshift=0;
scalar h_logloser=0.0398644;
scalar h_number=1.0986123;
scalar h_loglength=5.54;
scalar h_pceyr_1=64;
scalar h_pceyr_2=88;
scalar h_constant=0;

	generate x_betahat0 = SN_b1*h_logshift
							+ SN_b2*`a'
							+ SN_b3*h_logshift*`a'
							+ SN_b4*h_number
							+ SN_b5*h_loglength
							+ SN_b6*h_pceyr_1
							+ SN_b7*h_pceyr_2
							+ SN_b8*h_constant;
	generate x_betahat1 = SN_b1*(h_logshift+0.4534249)
							+ SN_b2*`a'
							+ SN_b3*(h_logshift+0.4534249)*`a'
							+ SN_b4*h_number
							+ SN_b5*h_loglength
							+ SN_b6*h_pceyr_1
							+ SN_b7*h_pceyr_2
							+ SN_b8*h_constant;
	
	gen prob0 =normal(x_betahat0);
	gen prob1 =normal(x_betahat1);
	gen diff=prob1-prob0;
	
	egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost2 (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
    
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
    
    local a=`a'+ 0.001;
    
    display "." _c;
    
} ;

display "";

postclose mypost2;

restore;

clear;

use pwsims2.dta;

gen MV = _n-1;
replace MV=MV/1000;

twoway (line diff_hat MV if MV<=0.28, sort) (line diff_lo MV if MV<=0.28, sort lpattern(dash)) (line diff_hi MV if MV<=0.28, sort lpattern(dash)), ytitle("Change in Pr(War), Min to Max Shift") yline(0, lcolor(gray)) xtitle("Losing Side's Total Military Capabilities") legend(off);
graph save Graph figure5a.gph, replace;

save pwsims2.dta, replace;
